# konrad.smolinski@gmail.com
# date: 26/11/2010
# last update: 29/11/2010
#
# info:
# 	DGP for binary outcome Y and discrete X={x_1,...,x_K}
# required:
# library: mvtnorm
# ------------------------------------------------------------------
rcDgp <- function(a,b,d,xthres,xval,zval){
	xthres <- c(-Inf,xthres,Inf)
	lambda <- d[2]^2 + 2*d[2]*d[3]*b[1] + (b[2] + b[1]^2)*d[3]^2 + d[4]^2

	K <- length(xval)
	R <- length(zval)
	rho <- list()

	for(r in 1:R){

		prY0Xk <- rep(0,K)
		prY1Xk <- rep(0,K)

		for(k in 1:K){
			mn <- c(0,0)
			covWQk <- d[2] + b[1]*(d[3] - d[2]*xval[k]) - d[3]*xval[k]*(b[2] + b[1]^2)
			varQk <- (1-xval[k]*b[1])^2 + b[2]*xval[k]^2
			sigMat <- cbind(c(lambda,covWQk),c(covWQk,varQk)) 

			prY0Xk[k] <- pmvnorm(lower=c(xthres[k]-d[1]*zval[r],-Inf),upper=c(xthres[k+1]-d[1]*zval[r],a[1]+a[2]*xval[k]),mean=mn,sigma=sigMat)[1]
			prY1Xk[k] <- pmvnorm(lower=c(xthres[k]-d[1]*zval[r],a[1]+a[2]*xval[k]),upper=c(xthres[k+1]-d[1]*zval[r],Inf),mean=mn,sigma=sigMat)[1]
		}

		# Pr(Y,X|Z=z[r])
		rho[[r]] <- rbind(prY0Xk,prY1Xk)
		colnames(rho[[r]]) <- paste("x",1:K,sep="")
		rownames(rho[[r]]) <- c("y0","y1")
	}
return(rho)
}

#Validate
# ------------------------------------------------------------------
valDgp <- FALSE

if(valDgp){
cat("\n Validate: rcDgp() : \n\n")
	a <- c(0,1)
	b <- c(0,1)
	d <- c(2,3^(-1/2),3^(-1/2),3^(-1/2))
	zval <- c(-1,-0.5,0.5,1)
	xthres <- c(0)
	xval <- c(0,1)

	dgp1 <- rcDgp(a,b,d,xthres,xval,zval)
	dgp1
}